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Abstract. We develop a formalism that can be used to model slowly rotating superfluid Newtonian neutron stars. 
A simple two-fluid model is used to describe the matter, where one fluid consists of the superfluid neutrons that 
are believed to exist in the inner crust and core of mature neutron stars, while the other fluid is a charge neutral 
conglomerate of the remaining constituents (crust nuclei, core superconducting protons, electrons, etc). We include 
the entrainment effect, which is a non-dissipative interaction between the two fluids whereby a momentum induced 
in one of the fluids will cause part of the mass of the other fluid to be carried along. The equations that describe 
rotational equilibria (i.e. axisymmetric and stationary configurations) are approximated using the slow-rotation 
approximation; an expansion in terms of the rotation rates of the two fiuids where only terms up to second-order 
are kept. Our formalism allows the neutrons to rotate at a rate different from that of the charged constituents. 
For a particular equation of state that is quadratic in the two mass-densities and relative velocities of the fluids, 
we find an analytic solution to the slow-rotation equations. This result provides an elegant generalisation to the 
two-fiuid problem of the standard expressions for the one-fluid polytrope £ oc p^. The model equation of state 
includes entrainment and is general enough to allow for realistic values for, say, mass and radius of the star. It 
also includes a mixed term in the mass densities that can be related to "symmetry energy" terms that appear in 
more realistic equations of state. We use the analytic solution to explore how relative rotation between the two 
fluids, the "symmetry energy" term, and entrainment affect the neutron star's local distribution of particles, as 
well as global quantities as the Kepler limit, ellipticity, and moments of inertia. 

Key words, neutron stars - rotation - superfluidity - entrainment 



1. Introduction 

Puls ars w ere first discovered over 34 years ago (Hewish 
et al 196^ ). Since then, nearly 1300 such rapidly rotating, 
highly magnetised neutrons stars have been found, and as 
pointed out by Lorimer in his recent Living Reviews in 
Relativity (Lorimer 2001), 700 of these were found in the 
last 4 years alone. Without doubt, the overwhelming ma- 
jority of these pulsars are old and cold (with core temper- 
atures below 10^ K). According to equation of state calcu- 
lations, their interiors will contain superfluid neutrons, su- 
perconducting protons, a plasma of highly degenerate and 
ultra-relativistic electrons, and perhaps other more exotic 
particles (pions, hyperons, etc.) deep in their cores. In re- 
cent years, there has been a continued effort to model su- 
perfluid dynamics in neutron stars in both the Newtonian 
(Epste in |1988| ; Mendell & Li ndblom |199lt Men deU |1991a 
1991bt Lin dblom & Mendell |l994 |1995|, |200C| ; Lee |1995 



This paper is aimed at improving our understanding of 
local and global properties of rotating Newtonian super- 
fluid neutron stars. 



The strongest evidence for superfluidity in the inner 
crust and core of a neutron star (see, for example, Sauls 
( 1989| ) and references therein) is provided by the well- 
known glitch phenomenon (the occasional sudden spin-up 
of relatively young pulsars). Our confidence in an expla- 
nation based on the transfer of angular momentum be- 
tween two loosely coupled components is bolstered by the 
fact that the neutrons and protons are described using the 
same many-body theory of Fermi liquids and BCS mech- 
anism that has been so successful at describing supercon- 
ductors (Pines 



Prix 1999; Se drakian & Wasserman ^000 ; Andersson & 



Com er 2001b ) and gene ral relativisti c regi mes (Carter 



1989|, C arter & Langl ois [19954 |l995b| , |1998[ La nglois e t 
al |1998| ; Comer et al [l999[ Andersson & Comer |2001aD . 



Nozieres 1966). This being the case, one 
would expect superfluidity in neutron stars to share many 
of the well-established properties of laboratory superfluids 
and superconductors. For instance, it is known that in a 
mixture of two interpenetrating fluids there is a coupling 
that arises whereby the momentum of one of the liquids 
is not simply proportional to that liquid's velocity, rather, 
it is a linear combination of the velocities of both fluids. 
This is the so-called entrainment effect and it implies that 
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when one liquid starts to flow then it will necessarily in- 
duce a monicntum in the other constituent. Entrainment 
between protons and neutrons is a key component of mod- 
els of neutron star superfluidity, and it is something that 
we will focus on in this study. 

We develop a formalism for describing slowly rotat- 
ing Newtonian superfluid neutron stars which allows the 
neutrons and protons to rotate at different rates. Our 
analysis is based on a Newtonian model that is the non- 
relativistic limit of a comprehensive model developed by 
Carter, Langlois and thei r coll aborators for the general 
relativistic r egime (Carter 1989 Carter & Langlois 1995a , 



1995U |1998[ Langlois et al |l998D . Our model is simplified 



in the sense that it describes a superfluid neutron star in 
terms of only two fl uids ( we refer the reader to, for ex- 
ample. Comer et al ( 1999|) for justification). One fluid is 
composed of the superfluid neutrons, existing in the inner 
crust and core, and the other fluid is a charge neutral con- 
glomerate of the remaining constituents (i.e. crust nuclei, 
core protons, and the crust and core electrons) that we 
will loosely refer to as the "protons." 

A further simplification of our model regards the vor- 
tices of the superfiuid. Because the density of neutron vor- 
tices is expected to be about 10^ — lO'^cm^^ for typical 
pulsar rotation rates, it is reasonable to adopt a smooth 
averaged description of vorticity on a macroscopic scale. 
The resultant model then resembles a perfect fiuid. There 
is, however, one important difference that arises because 
of the possible interaction between the vortex lattice and 
the second fiuid. This is a dissipative effect usually re- 
ferred to as "mutual friction" . In general, this mechanism 
tends to drive the two fluids towards co-rotation, but there 
are two extreme cases where a stationary description of a 
two-fluid star with different rotation rates would still be 
possible. The first case is (obviously) when the interaction 
is very small (the "free" vortex limit). Somewhat surpris- 
ingly, similar results apply in the case when the interaction 
is very strong (resembling the "pinned" vortex case in a 
solid crust, cf. Langlois et al. 199^ ). Although we could, 
in principle, allow for this second case in the present con- 
text (as has been done in an earlier study by Prix 1999| ), 
we will refrain from doing so and assume that vortex fric- 
tion is negligible on the time-scales we are considering. 
This is done because the emphasis of the present work is 
on the detailed role of the entrainment effect. Including 
vortex friction would unnecessarily complicate the discus- 
sion and distract the attention away from the new piece 
of physics we have incorporated. 

We will generalize the Chandrasekhar-Milne slow- 
rotation method (Chandrasekhar 1933| ; Milne 1923) to our 
simple two-fluid model. That is, we make an expansion in 
terms of the fluid rotation rates keeping only terms up 
to second order. We include entrainment, and consider 
the most general case where the neutrons do not coro- 
tate with the protons. This relative rotation of neutrons to 
protons is also present in the general relativ istic slow rota- 
tion scheme of Andersson & Comer ( ^OOla ) (who provide 
an extended discussion as to why such a generalization 



is necessary), and the Newtonian slow rotation scheme of 
Prix (199£). However, although the formalism developed 
by Andersson & Comer does allow for entrainment, their 
numerical study did not include it, and Prix excluded en- 
trainment altogether. As we wish to explore the effects 
of entrainment on rotational equilibria, and the current 
best model for entrainment is for the Newtonian regime 
(Borumand et al 1996), it is natural for us to take up where 
Prix (1999) left off, and consider Newtonian models. 

As an application of the slow-rotation formalism, we 
will consider a particular form of the equation of state 
(EOS) that is the most general quadratic form in the mass- 
densities and in the relative velocity. This EOS contains a 
mixed term in the neutron and proton mass densities that 
can be related to so-called "symmetry energy" terms in re- 
alistic equations of state (i.e. terms that vanish when there 
are equal numbers of protons and neutrons). Remarkably, 
we find an analytic solution to the slow-rotation equations 
for this EOS. This solution has enough free parameters 
that realistic neutron star values for the mass, say, can 
be obtained. Using this solution we will show how rela- 
tive rotation, entrainment, and the "symmetry energy" 
term affect the neutron star's local distribution of par- 
ticles, as well as global quantities like the ellipticity, the 
Kepler limit, and the moments of inertia. 

Before concluding this introduction we should note one 
of the main motivations behind the present work. Once 
we have developed a framework for obtaining stationary 
superfluid models, we want to study their dynamics as 
manifested by the various modes of pulsation. This issue 
is of particular interest as it is known that several different 
modes of a rotating neutron star are generically unstable 
due to the emis sion of gravitational radi ation (Friedman 
& Schutz 1978 ; Andersson & Kokkotas 2001 ). One can 
even plausibly speculate that pulsation modes unique to 
superfluids will be excited during a pulsar glitch and that 
the resultant gravitational radiation may be detected by 
a future generation of advanced detectors (Andersson & 
Comer 2001c ). Our aim is to investigate these various 
possibilities by using the models developed here as back- 
ground for a linear perturbation calculation of the relevant 
pulsation modes of a rotating superfluid neutron star. 

2. Canonical two-fluid hydrodynamics 

In this section we present the Newtonian canonical de- 
scription of a non-dissipative interacting two-fluid system, 
which can be derived either as the non-relativistic limit 
(Andersson & Comer 2001b| ) of the canonical covariant 
description developed by Carter and Langlois (1998), or 
directly from an analogous Newtonian variational princi- 
ple (Prix 2000, 2001). Since we are mainly interested in 
describing a superfluid neutron star core we label the two 
fluids by indices "n" and "p" , representing the neutrons 
and protons, respectively. The fundamental variables of 
our description are then the respective particle number 
densities «„ and np and the corresponding particle cur- 
rents rin and rip. We will only consider situations where 
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the particle numbers are conserved individually, so that 
we have 



dtUn + Viul = , and dtUp + V; 



0, 



(1) 



where we sum over repeated spatial indices i,j = 1,2,3. 
The respective velocities t>n, I'p and the mass densities p^ 
and pp follow from 

(2) 



n„ = 


nnVn , 


Up 


= npVp , 


and 








Pn = 


m"nn , 


Pp 


= rTi^rip 



P = Pn + Pp: 



(3) 



where m" and m^ are the respective (fixed) masses per 
particle of the two fluids, and p is the total mass density. 
The internal energy density or "equation of state" 
£{nn,np,nn,np) of the two- fluid system must satisfy 
Galilean invariance (and isotropy, as we want to describe 
isotropic fluids), and therefore has to be of the form 



£ = £{nn,np,A'^) , where 



A2 



^p)' 



(4) 



This energy function defines the respective chemical po- 
tentials /^", pP and the "entrainment function" a as the 
conjugate variables to Um rip and A^, namely 



d£ — p^ dn-a -\- p^ dup 



.dA^ 



(5) 



which can be regarded as the "first law of thermodynam- 
ics" for this system. The canonical description of the two- 
fluid system is based on a convective variational princi- 
ple (Prix 200C,2001) analogous to the method used by 
Carter () in the general relativistic context. The dynam- 
ics of the system is therefore described by a Lagrangian 
density >C(nn,rip, rini'^^p), which has the usual form of 
"£ — kinetic energy — potential energy" . In the present 
context this means that 



^ = l^Pr^vl 



Pp^p - (^ 



P*), 



(6) 



where $ is the gravitational potential, which is related to 
the total mass density p by the Poisson equation 



V2$ = 47rGp. 



(7) 



Variation of the Lagrangian density C defines the mo- 
menta per (fluid) particle, p" and pP, and the "rest frame 
chemical potentials," — pg and — Pg, as the conjugate vari- 
ables to the currents rin, rip and particle densities n„. 



'p? 



and the conjugate momenta p^ and pP are related to the 
particle currents as 



p" = /C™nn+/C"Pnp 
pP = /C"P nn + /CPP np , 



(10) 



where the symmetric "entrainment matrix" /C has the 
following components: 

/C°" = \ (m"n„ - 2a) , /CPp = ^ {ivFup - 2a) , 



/C"P = 



2a 



(11) 



The general form of the entrainment relation (|lO[), 
which states that the momenta are linear combinations of 
the currents, is a consequence of the Galilean invariance 
of the internal energy £. The most important difference 
from the case of a single fluid is that the particle momen- 
tum of each of the two fluids is in general not aligned with 
its respective current. From (10) we see that this devia- 
tion is caused by the off-diagonal matrix element /C"p of 
(0), which is proportional to the entrainment function 
a and therefore expresses the dependence of the internal 
energy density £ on the relative velocity A, as seen in 
(ra). Intuitively this means that a particle current n^ of 
the neutrons impinges some momentum on the protons 
and vice versa, due to the interaction between the two 
fluids. In the absence of such an interaction, i.e. for a = 0, 
the general entrainment relation (IlG) simply reduces to 
the usual single-fluid form, and we have p" = m^Vi^ and 



P' 



mPvp. 



The equations of motion for the two fluids, obtained 
from a Newtonian convective variational principle (see 
Prix pOOOl pOOll Andersson & Comer pOOlbj ), thus have 



the following form: 



n„(9ip" 
"-p(9tpP 



Vpg) 



^i{V,p^~Vp^) 



(V,pP 



VpP) = 



(12) 



I.e. 



At this point, the two equations are "formally" uncou- 
pled, but the interaction between the two fluids enters 
through relations (g) and (|lO|), which relate the kinemat- 
ical quantities n^, rip, rin and rip to their respective con- 
jugate dynamical quantities pg, pg, p" and pP. In fact, 
the fluids will be coupled even in the case of two "non- 
interacting" fluids, i.e. for an equation of state of the form 
£ = £n{nn) + £p{np). This is simply because the t wo flu ids 
live in the same gravitational potential, cf. Prix (1999). 



dC — p^ ■ dn,,i + pP • dUp + Pgrfrin + Pgdrip . 



Using the explicit expression (g) for the Lagrangian den- 
sity together with the "first law" (g), we obtain the fol- 
lowing expressions for Pq and Pq , 



Po = p" + v,,-p' 



-Pi- 



Vp-p' 






m"$ , 



mP$, 



^ ^ 3. Entrainment and effective masses 

In the framework of condensed matter physics the descrip- 
tion of interacting independent constituents, for example 
electrons moving through a "background" of ions, is usu- 
ally formulate d in t erms of an "effective" or "dynamical" 
mass (Ziman 1965 ), rather than in the entrainment for- 
malism presented in the previous section. One reason for 
(9) this might be that in the usual contexts one of the two 
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constituents (the "background") can usually be assumed 
to be at rest with respect to the observer. It is then conve- 
nient to formulate the dynamics of the second constituent 
in a form that resembles the vacuum expressions, with the 
free particle mass being replaced by an "effective mass" 
that accounts for the interaction with the background. 

The general definition of the effective mass is based on 
the particle energy spectrum, E{p), say, of particles with 
momentum p moving relative to the non-vacuum back- 
ground. The velocity v of the particles as a function of 
their momentum is then given by t; = dE/dp, and the 
acceleration i) (where the dot represents the total time 
derivative) is linked to the force p via a relation of the 
form 



Pj 



where 



dpidpj 



(13) 



In the two- 



which defines the effective mass tensor i 
fluid problem we are interested in, the background can 
usually be assumed to be isotropic, and therefore the par- 
ticle energy spectrum will be of the form E = E{p'^), which 
implies that 



^dE 
dp"^ 



P- 



(14) 



In this case the effective mass tensor of ( [I4 ) can be ex- 
pressed as 



9^^'lA 
' V ' '' 



^ d^E 

dp^dp^ 



PiP] 



(15) 



In the limit of small momentum, where the energy will be 
mainly quadratic in p, we can neglect the second term and 
define an effective mass scalar m* in the usual way. This 
leads to 



const . 



(16) 



— = 2j— ^=> E{p) « -^ 
TO* op'' 2m* 

We see from (O) that an equivalent form of this defi- 
nition is 



P 



(17) 



which clarifies the link between the effective mass and the 
more general formalism of entrainment used in the present 
work. The key observation to make is that the effective 
mass description is formulated in the "rest-frame" of the 
background. If we choose the background to be the "neu- 
trons," we could define the neutron rest-frame by rin — 0, 
and therefore the entrainment relation ( |lO| ) in this frame 
reads as 



P' 



(/CPPnp) v^ 



(18) 



By comparison with (llj) we can relate the proton effective 
mass mP* to the entrainment matrix by mP* — KP^n-p. 
Using the explicit expression dnl) for /Cpp in terms of the 
entrainment function a, we find 



2a — np{imP — nv^*) — epp . 



(19) 



where we introduced the dimensionless "entrainment co- 
efficient" e asR 

mP — mP* 

e ^ . (20) 

mP 

Therefore the entrainment matrix /C^^ is determined com- 
pletely in terms of the effective mass of a single particle 
species, for example the protons in the present formula- 
tion. Of course, we also see (from (19) by symmetry) that 
the effective masses of the two species are not indepen- 
dent, because they obviously have to satisfy 



m 



^(toP 



^)- 



(21) 



However, as pointed out by Carter ( pOOl , private com- 
munication) , the above definition of the effective mass is 
not quite unique because we have to specify what we mean 
by the "rest-frame" of the neutrons. Indeed, we have de- 
fined it by setting Tin = 0. However, we see from (ll^) that 
in this frame we generally have p" ^ 0. Therefore another 
equally viable choice of the "neutron rest-frame" would be 
given by setting p" = 0. This would then lead to nn 7^ 0. 
Based on this second choice, the entrainment relation (10) 
becomes 



pP= /cpp- 



np^2 



(/C"P ) 



(22) 



leading to an alternative definition of the effective proton 
mass mP'^ , which is not equivalent to that of Eq. (ttq). The 
corresponding relation between a and toP"*^ is 



2a = np{'mP 



.P*\ 



) 1 + 



Up (mP 



.P# 



m" 



(23) 



How does this apparent ambiguity in the definition of 
the effective mass affect our attempt to model entrainment 
in superfluid neutron stars? Fortunately, it turns out that 
the difference between the two (equally valid) definitions 
becomes numerically small in two interesting special cases: 

(i) neutron star matter: toP ^ m^'^ ~ tti" and rip ^ rin, 



(ii) electrons in a metal: Ji^ = n+ and rri- 



m_ ^ m^ 



where in the second case n_ and n+ represent the number 
densities of negative and positive charges respectively (i.e. 
electrons and protons), and to_ and to_|_ the correspond- 
ing mass per particle. In these two cases equation (23) 
becomes approximately equal to the previous expression 
(|l9|), and the two possible definitions lead to very similar 
results. 

There exist in the nuclear physics literature a few cal- 
culations of the proton effective mass at neutron star den- 
sities. The results are very much equation of state de- 
pendent, but one can extract some useful constraints on 



We note that Lindblom & Mendell (200C) have used a 



slightly different parameter to characterize the entrainment ef- 
fect. Their parameter e is related to our coefficient e via 



ep 
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the range of values. Chao et al (1972) find, for the range 
pnuc ^ P ^ 2.5/9nuc (where pnuc represents nuclear density, 
i.e. Pnuc ^ 2.7 X 10^'^ g/cm^), that m*Jmp « 0.6 - 0.5; 
Sjoberg ( |1976| ) finds, for pnuc < p , that m^/mp « 0.6-0.4; 



and more recently, Baldo et al (1992) find, for the same 
range as Chao et al, that m*/m-p k, 0.7. Thus, we see that 
the proton effective mass apparently can range over the 
values 0.3 < m*/mp < 0.7. The calculations also show 
that the effective proton mass varies slowly with the den- 
sity. This means that one can, as a rough but reasonable 
approximation, assume that the effective proton mass re- 
mains constant throughout the core of a neutron star. This 
is, in fact, the case for the analytic solution that will be 
discussed later. 

4. Stationary two-fluid neutron star models 

Stationary two-fluid configurations are characterised by 
9tp" — and dtp'^ = 0, and therefore the equations of 
motion (|T2|) take the form 









rrf^v^ 



mPv^ 



(p^V,«i + <V,p^) = 0, 
(pPV.z;^-F«^V,pr) = 0(24) 



where we have inserted the explicit expressions (g) for pg 
and Pq. We consider stationary configurations where both 
fluids are rotating uniformly with rotation rates Hn and 
rjp about a common axis (otherwise the system would not 
be stationary). 



< = f7„(^' , and vl = npif' , 



where tp'^ is the assumed axial symmetry vector with norm 
{'fi^ipi)^!'^ = w, and w is the distance from the rotation 
axis. In spherical coordinates (r, 9, (p) we, of course, have 
w = rsvaO. As both fluids are flowing in the same direc- 
tion (/?*, we can write 



^ ^P ^ 



(26) 



which allows us to express the momenta p" and pP of (|10| 
simply as 



P' 



P' 



M"t;, 



with M" 



f7r 



MP^n, with MP = /CP"nn+/CPPnp--ii. (27) 

i in 



Because of the axial symmetry, scalars like Af" and 
Af P cannot depend on the angle (p, and therefore we have 



viVjAP' = , and vlVjAP = . 



(28) 



Furthermore, since we are assuming uniform rotations, cf. 
(Eq), we can deduce the identity 



<VjW„i 



1 



v.^^' 



^nl 



(29) 



This means that the second parenthesis in the equations 
of motion ^ vanishes. Namely, using (|2|), (||) and (H), 
we find 

p^"V,«^ + vrsjjp- = Af " Qv,«2 ^ „^^v,z;n.) = , 

This reduces equations (p3) to two first integrals of mo- 
tionn, namely 



p" + m"$ m"fi2 ^ Qn ^ 

2 

pP + mP$-^mP02 ^ ^P^ 



(30) 



where the constants C^ and C^ are in general determined 
by the rotation rates U,^ and fip. These first integrals, to- 
gether with Poisson's equation (R), an equation of state 
£{n-m Tip, A^) and appropriate boundary conditions com- 
pletely specify the solution. Note that these first integrals 
are also the Newtonian limits of the gene ral rela tivistic 
results obtained by Andersson and Comer (2001a). 

It will be convenient in the following to write our 
equations in dimensionless form by expressing all quan- 
tities in their "natural units." We choose these to be 
the radius R of the (non-rotating) star for lengthscales, 
its central density pQ for densities, and 1/\/AttGpq for 
timescales. For simplicity of notation we will continue to 
use the same symbols for the dimensionless quantities, so 
Poisson's equation (|7|) now reads 



(25) V2$(r,0)-p(r, ^) 



(31) 



whereas the form of the flrst integrals (pffl) is unchanged. 
We note that the "natural rotation rate" \/47rGpo is 
slightly greater than the maximal rotation rate, the Kepler 
limit, f^K at which mass shedding will occur. 

For realistic equations of state the mass shedding 
limit is well approximated by the empirical expression 
ilK ~ \/47rGp/3, where p is the mean density. Therefore 
the dimensionless rotation rates fin and fJp are generally 
smaller than 1, and for typical rotation rates fi of most 
observed pulsars we will have 



f^<l. 



(32) 



Note that the difference between the "natural rotation 
rate" and the mass shedding limit is further amplified 
by the fact that po ~ 5 — lOp for a typical neutron star 
equation of state. That being said, we must still proceed 
somewhat cautiously near the Kepler limit. A compar- 
ison between slow-rotation results and numerical calcu- 
lations for rapidly rotating stars show that, while the 
slow-rotation approximation can accurately describe the 
fastest observed pulsars it deteriorates significantly near 
the Kepler limit. This is illustrated in Fig. |l|, where we 

^ An elegant and more generally applicable method of ob- 
taining these first integrals is presented in Appendix \p\. 
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Fig. 1. Comparison of one-fluid slow-rotation configura- 
tions with numerically determined one-fluid configurations 
(obtained using the LORENE code). It is clear that the 
slow-rotation approximation leads to a significant under- 
estimate of the rotational flattening for stars rotating near 
the Kepler limit. 



compare the slow-rotation results for the equatorial and 
polar radii of ordinary, one-fluid stars (governed by an 
equation of state of the form £ (x p^) to accurate numerical 
results. The latter were obtained using the LORENE code, 
based on pseudo-spectral techniques, which was developed 
by the M eudon Numerical R elativ ity group (Bonazzola et 
1993 ; Gourgoulhon et al. 1999| ) and which can be used 



al 

to build rapidly rotating, Newtonian and general relativis- 
tic neutron stars. We see that the slow-rotation approx- 
imation works quite well up to and including the fastest 
known pulsar, but starts to fail (by some 15% to 20%) as 
the Kepler limit is approached. These arguments motivate 
the range of applicability of the "slow-rotation expansion," 
which we will use in the following sections. 



5. The slow-rotation approximation 

5.1. Slow-rotation expansion of the two-fluid model 



In order to approximate the solution to (pOD and (|3l| 



we apply a method initially due to Chandrasekhar ( [1933D 
and Milne ( 1923| ). This method is based on the assumption 
that the rotation is slow enough that the configuration is 
only slightly oblate, in such a way that one can express it 
in terms of a Taylor expansion around the non-rotating 
spherical solution. In the following we consider such a 
"slowly-rotating" two-fluid star. In the previous section 
we have seen that this corresponds to f2n and fip being 
small compared to the natural rotation rate ^/A^^Gpo. In 
the dimensionless form of the equations, the small param- 
eters of the expansion are simply the rotation rates, as 
J7n '^ 1 and lip ^ 1. We can therefore write the solution 
as a Taylor expansion in il^ and flp. In doing this we will 
neglect all terms beyond second order. 



Under these assumptions any scalar physical quantity 
Q of the rotating star can be written 






dQ 



1 d^Q 

2 dni 

1 d^Q 



an. 



fir 



2 dni 



ni 



ni 





OH 






U Lr f 1 U U "p 



(33) 



Furthermore, because we are actually expanding in the 
"vector" (f2n,r2p), it will be convenient to introduce the 
notation Q for the "constituent vector" with components 
fix, where X = n, p is the constituent index. Thus we have 



n = {fin, ftp) , i.e. (H) 



'X 



ft 



■'X 



(34) 



Terms of odd-order in fl in the expansion (p3) will vanish 
identically because the configuration has to be invariant 
under a simultaneous inversion of both rotation rates, i.e. 
for fin — > —fin and flp —^ —flp. 

We know that the static solution Q\n=o is purely spher- 
ical and denote it Q*-°-'(r'), i.e. 



Q{r,9;n = Q) 

With these definitions, the second-order expansion 
Q can be written in a more compact form as 

n)^Q'^"\r)+SQir,0)+O{ff) 



Qir, 
with 
SQ = flxQ^^fl 



(35) 
3l)of 

(36) 

(37) 



where we automatically sum over repeated constituent 
indices (this convention is assumed from now on), and 
Q^^ (r, 9) is the symmetric expansion coefficient matrix 



Q^'^{r,0) 



d^Q 



dflxdfl^ 



(38) 



n=o 



It has been stated severa l times in t he literature 
(Mon aghan & Roxb urgh |1965| ; Hartle |1967| ; Smith 



1975 



1976 ; Tassoul 1978 ) that the slow-rotation expansion of 
the density p{r, 9) breaks down when approaching the sur- 
face of the non-rotating star, i.e. for r —^ 1, because there 
p(o)(l) = 0. However, it is evident from ( p6|) that the va- 
lidity of this expansion is completely independent of the 
values oiQ'-'^\r). The only condition for its validity is that 
the neglected terms 0{fl ) are effectively small compared 
to the 0{fl^) terms included in the analysis. 

We now expand all physical quantities of (pfl) and (31 ) 
up to second-order in fl: 



<P{r,9) 


= $("' +5$, 


(5$ = fix <^^ fly , 


(39) 


p{r, 9) 


= p(°)+<5p. 


6p = nxf^''nY, 


(40) 


nA{r,9) 


= n\' + Sua , 


Sua = fix nY ^y , 


(41) 


p\t,9) 


= /(°) + 5/ , 


(5/ = fix /'^^ fly , 


(42) 


C^ 


= c^(") + 5C^ , 


5C^^flxC^'''''flY. 


(43) 
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The total mass density p is, of course, linked to the indi- for a given equation of state £(n^, n-p, A^). Thus we define 
vidual particle densities by (H), i.e. p = m^ ua- Therefore the "density structure function" Sxrir) as 
we have the relations 

' . d^£ 
and ff^^ 



p(0)=m^4°) 



m^nY 



(44) ^^"^^'^={^1) 



dnxdny 



(55) 



Inserting the various expansions into ^ and (|l|), the ^^ich is symmetric (and we assume invertible with in- 



zeroth-order equations, which determine the non- rotating 
configuration, are found to be 



\XYs 



verse {S ^) ), and the "entrainment structure func- 
tion" ^{r) as 






while the second-order coefficients of 

,2 



(45) 
(46) 



/3^(r 



dp^ 



dA^ 



d^£ 



dnxdlS? 



da 



dnx 



(56) 



satisfy 



M 



m 



T 



A.XY , „,A^Y _ ^ rA,Xy _ ^A,XY 

A^XY 



V^^^ = 



m 71^ 



(47) 
(48) 



where we used the definition (^) of the chemical potentials 
p^ and the entrainment function a. Comparing (5^) to the 
expansion ( |42| ) we can identify 



/(°) =/(n(°),n(°),0). 



In equation (B7h we have introduced the constant matrices A* 
5^'-^^ , which are defined as 



,.A.XY 



= S 



-l\^ „XY 



zu^p^A 



XY 



(57) 
(58) 



srn,XY 



to" 




6P 



,XY _ 





TOP 



and we have thus arrived at an algebraic relation between 
the p^'^^ and the r?^^. 
(4^) Inserting this relation into the &-st integral (p7|), we ob- 

tain an explicit expression for the density coefficients n^^ 



These allow us to write the second-order terms rn^Vll^ and in terms of the gravitational potential coefficients $^^ 



mPfip m 



OD as 



Specifically, we have 



and 



^x 5'''^^ ny 



mPnl = nx 6^'^^ fiy 



5.2. Reduction to a single equation 



(50) 
(51) 



■n^^ 

n^ 



5. 



'AB 



a 



i,XY 



^(^^ 



,XY 



,B aXY\ 



2(3" A- 



^B^rn 



(59) 



We can find an algebraic relation between the chemical po- 
tential coefficients p/^'-^^ and the density coefficients n^^, 
which will allow us to considerably simplify the problem. 
This relation is found by expanding the chemical poten- 
tial p/^{nn,np,A^) defined in (||) in its arguments up to 
second-order in fl. Using ( |4l| ) we find 

/(nn,np, A2) = /(40),n("),0) 



Introducing the "derived" background functions 

E^'-ir)^ i5^(r)(5^^^^-2/3^(r)A^^), 

kA(r) = SAB{r)m^ , 
equation ( |5S| ) can be written 






(60) 



(61) 



O. 



/ V 



\dnz 



„XY 



which is a reformulation of the first integrals of motion 
( p7| ) in terms of the second-order coefficients. 

The total mass density coefficient /9^^ can now be writ- 
ten 



9/ 



9A2 



^2^XY 



Uy + 0{n^) , (52) 



/^ 



m^nY 



3n7 



where we have defined the constant matrix A''''^ as 



= TO^6kB(r)C^'^^ + :^ii;^^(r)-fc(r)$^^, (62) 



vxy 



1 -1 
-1 1 



where we have introduced the further abbreviations 



(53) 



I^^{r) ^ m^E^^'ir 



This definition allows us to write the relative velocity 
squared, namely A^ = ti7^ (fl^ ~ ^p) i i^i the form 



1 



w'^ nx A^^ ny 



(54) 



- -,«^5^(r)(5^^^^-2/3^(r)A^^) 
k{r) = TO kA{r) — m SAB{r)m , 



(63) 



which are "known" functions of the non-rotating star's 

The partial derivatives in (|5j) are evaluated for the configuration, determined in terms of the "structure func- 

non- rotating configuration, and thus they are assumed to tions" iSxy(f) and (3^{r), together with the constants 

be known functions, depending only on the static solution S"^'^^ and A^^ and the constants of integration C"^'-^^ . 
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These constants are determined by the boundary condi- 
tions. Inserting (p2|) into Poisson's equation (B^ reduces 
the problem to a single partial differential equation for the 
coefficients ^^ (r, 6) , namely 



3n7 



V2<i>^^ + fc(r)*^^ - :^E^y{r) + m''SAB{r)C'^ 



,XY 



(64) 



Given the non-rotating background, a solution to this 
equation fully specifies a slowly-rotating, two-fluid con- 
figuration. 

We note that a necessary condition for this method to 
work, i.e. for ( |6^ ) to be well-defined, is that the struc- 
ture function SxYif) is regular everywhere inside the 
star. However, it is well-known that in the case of a 
single fiuid the standard polytropic equations of state, 
i.e. £ = Kn^+^/^, leads either to a vanishing (for iV > 1) 
or an infinite density gradient (for iV < 1) at the surface. 
This behaviour could make ( |55| ) singular at the surface. In 
order to clarify this point, we can express the definition 
( pq ) equivalently as 



.f'(r)=5.>-M^(°)'(r), 



(65) 



where a prime denotes the radial derivative d/dr. It is 
seen from ([45| ) that the chemical potentials ^'^"•'(r) and 
their derivatives are regular everywhere, even at the sur- 
face. Therefore we can conclude that the structure func- 
tion SxY behaves like the density gradients n^ at the 
surface. The slow rotation expansion used in the present 
work is therefore well-defined for stars with finite or zero 
density gradients at the surface, but is not applicable to 
stars with an infinite density gradient at the surface, sim- 
ilar to the iV < 1 case for the polytropes, as the structure 
function Sxv then diverges at the surface. 

Inserting ( [45|) into (pq), we obtain expressions for the 
derived structure functions /c4(r) ( |60[ ) and k{r), which in- 
dicate their physical meaning. Using ( pOJ ) we have 



,(0)' 



-fcA$("^' 



dn 



(0) 



d$(o) 



and from the definition mdi) it follows that 



^(0)' 



-fc$(")' 






(66) 



(67) 



From these relations we see that the "structure" functions 
/cA(r) reflect the change of the individual number densi- 
ties with respect to the gravitational potential, while the 
function k{r) reflects the corresponding change in the to- 
tal mass density. These relations proved very useful in 
deriving the results discussed in Section 7. 



the orthogonal basis of Legendre polynomials Pi (6) . These 
are the eigenf unctions of the operator Vg, satisfying 



V^Pi{0) 



1(1 + 1] 



PiiO)- 



Therefore, writing 

oo 

leads to 
V2$^^=Vp,(0)A$f^(r), 



1=0 



E 

l=Q 



where the differential operator T>i is defined as 
d^ 2d l{l + l) 



V, 



dr^ 



■ dr 



(68) 



(69) 



(70) 



(71) 



which represents the radial part of the Laplacian with 
an additional "angular momentum" part proportional to 
1(1 + 1). 

Using ([70|) together with the identity 



3tn2 



^r'[l-P,{d)], 



(72) 



we arrive at a series of ordinary differential equations for 
the coefficients ^^(r) in (|69|). We have 



2 77-xy 



r^E 



m^'SABC' 



Vi^^+k^^^O, for />4 



■,XY 



;^5 



(73) 
(74) 
(75) 



where only terms of even I will be non-vanishing be- 
cause of the expected equatorial symmetry of the solution. 
Furthermore, we will see in the following section that the 
boundary conditions are homogeneous^ and therefore the 
homogeneous equations ([75[) have only trivial solutions. 



l>4: ^^{r 



0. 



(76) 



The only non-vanishing contributions to ^^ (r, 6) there- 
fore consist of the two functions $^^ (r) and $^^ (r) , and 
the problem has been reduced to that of solving the two 
differential matrix equations (|7^) and (JTJ). This is, of 
course, in complete analogy with the standard result of 
the single fiuid slow-rotation approximation. 

Once the solutions <I*f ^ to the differential equations 
[3[) and (^ have been obtained for a given equation of 
state (subject to the boundary conditions to be discussed 
in the next section) , the individual density coefficients n^ 
are determined from (m]). 



5.3. Separation of variables g Boundary conditions and integration constants 

So far, the slow-rotation approximation has allowed us ^ , o j j-j.- 

, 1 • • • 1 ,1 /CTTi , /TTHx • 1 0.1. Boundary conditions 

to reduce the mitial problem ( |3(J|) and (|31|) to a smgle 

equation ( |64|) for a single unknown function $^^(r, 9). We As Eqs. (JTg) are linear, second-order differential equations, 

now separate the variables r and 9 by expanding ^^ in two boundary conditions are necessary for each equation. 
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In the case of the I — equation (O) we need two further 
conditions in order to fix the two constants of integra- 
tion C"'"^^ and CP'^^. First, we require that the solution 
must be regular at the centre of the star. At r = the dif- 
ferential operator Vi is singular, which leads to the first 
boundary conditions 



and of the radial derivatives: 



/ > 2 



$f^'(0) 



0, 



<i>f^(0) = 0, and ci,f^'(0) = 0. 



(77) 



where the primes represent the radial derivative d/dr. 
Secondly, we require continuity of the gravitational po- 
tential and its derivative across the star's surface. That 
is, we have to match the interior solution <i>(r, 9) to the 
exterior gravitational potential ip{r, 9) at the surface of 
the rotating star R{9): 

^\R(e) - V-li^w ' and $'|^(g) = ^'1^(9) • (78) 

Here, the radial derivative deviates from the normal 
derivative with respect to the actual surface R{9) only 
by terms 0(11''), which have been neglected. 

We write the slow-rotation expansion of the exterior 
potential ^ in the usual way as 

V'(")(r) + <5V'(r, 9)+Oiif), 



with 

5V = ^x V^ny 



(79) 



(80) 



and expand the second-order coefficient in Legendre poly- 
nomials: 



^^(r,0) = ^V^^(r)P;(0), 



(81) 



1=0 



Because ip is the solution to the Laplace equation V^V — 
(normalised by lim^^oo V' = 0), we know the radial eigen- 
functions for all /, namely 



^^ir) 



^■XY 



rl+l 



(82) 



where the k^^ are constant coefficients. The continuity 
condition (kSh for the static solution implies 

$(o)(l)== V(0)(1) and $(")'(l)==:i/;(*')'(l), (83) 

and because the static potential <i>*^°^ at the surface also 
satisfies the Laplace equation, i.e. V^$^°^| =0, the 
second derivatives must also be identical: 

$(o)"(l)=^(")"(l). (84) 

Up to second-order in £2 the surface of the rotating star 
is 

R{9) = 1 + 6R{9) + O(il^) , (85) 

where we have used i?'"' = 1. This allows us to write 
the second-order expansion of the internal and external 
potentials at the matching surface R{9) as 

$|^(,) = $(o)(i) + ci>(o)'(i)^i? + <5$(l,0), (86) 

V'Ikw = ^*°Hl) + ^^°''(1) 6R + S^Pil, 9) , (87) 



$'|^^g^ = ^(''y(l) + <^(''y'{l)5R + 5<^'{l, 



(o)Vn j_ <ft(o)"('i^ AR J- Ads'n ft\ ^ (88) 

V^'lflce) = ^'^''^'(1) + ^^°>"(1) 5R + <5^'(1, 9) . (89) 



Therefore the matching conditions ( |78[) together with (|83| ) 
and (jsj) simply reduce to 

(5$(1, 61) = 5*(1, 6*) , and 5^'{\, 9) = 5^'{\, 9) , (90) 

which shows that the resulting boundary condition is — 
since there is no terms involving 5R — independent of 
the actual second-order surface of matching. It is thus not 
necessary to determine the perturbed surface of the star 
expficitly (by the condition /9|^/g^ = 0). Using (|8l|), (82) 
and ( |69| ) for the second-order terms 5ijj and ^$ in (90) 
directly results in the second boundary condition 



$f^'(l) + (? + l)$f^(l) = 0, for />0. 



6.2. Fixing the integration constants 



(91) 



As stated earlier, the I — equation requires two ad- 
ditional conditions in order to fix the constants of inte- 
gration C'^''^^ and CP'"^^. These conditions are crucial as 
they determine which type of rotating star sequence (as 
a function of the rotation rates iln and fip) the solution 
will describe. Among the many possible sequences, two 
seem particularly useful and will be described here. The 
first is the fixed central density sequence (FCD), which 
is probably the most straightforward and therefore most 
commonly considered choice. 

The FCD sequence is characterized by the simple con- 
dition (for ^ = n, p) 



nAL=n 



,(") 



r=0 



^4o(0) = 0. 



(92) 



This condition obviously implies that the respective total 
masses will change as functions of the rotation rates. This 
sequence therefore does not describe the same physical 
star at different rotation rates, but it has the advantage 
of leading to a very simple condition. From (p3) we find 
that we should require 



C^^^^=m^<I>^^(0), 



(93) 



The FCD sequence is, however, not particularly rele- 
vant from a physical point of view. For example, if we con- 
sider an isolated neutron star that is spinning down due 
to a magnetic torque we would expect the central density 
to increase as the star becomes less oblate. The physi- 
cally most interesting sequence of rotating stars thus cor- 
responds to requiring a fixed mass (FM). In other words, 
one has to impose the condition of constant respective to- 
tal masses (or equivalently, total particle numbers) of the 
two fluids if the sequence is to describe the same physi- 
cal star at different rotation rates. This means that (for 
^ = n, p) 



HAir, 9)dV 



n^^\r)dV. 



(94) 
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To second-order in the rotation rates, this condition re- 
duces to 



r nXf^^ir) dr = . 



(95) 



Inserting the explicit expression ( plj ) for the n^ q compo- 
nent yields the integral condition (for^,_B = n, p) 



[SabC 



.XY 



E^ - kA<^'^ (r)] dr = , (96) 



which allows us to determine the constants C^'^"^ for the 
fixed mass sequence. However, in practice one only needs 
to evaluate this integral for one of the two fluids, because 
the second fixed mass condition can be replaced by the 
equivalent but simpler requirement of fixed total mass, 
i.e. the integral over p — m^ nx, which analogously to ( |95[) 
leads to 



r^rrTn^^Q dr ^ . 
The / = component of Poisson's equation 

allows us to reduce this condition to 
$^^'(1)^0, ^ <i>^^(l) = 0. 



, I.e. 



(98) 



(99) 



where the second (equivalent) expression has been ob- 
tained from the second boundary condition (pi|). 

7. Ellipticities, moments of inertia, and Kepler 
rotation rate 

Before proceeding to discuss results for a specific model 
equation of state it is useful to digress on what kind of in- 
formation we want to extract from the calculation. There 
are obviously many alternative ways of describing rotating 
configurations. We will focus our attention on the elliptic- 
ities, the moments of inertia and the Kepler rotation rate. 
The first and second are interesting because they highlight 
how rotation affects the shape and mass-distribution of 
the star, while the last describes the limit of rotation at 
which mass-shedding at the equator sets in. 

We can explicitly express the respective fluid surfaces 
Ra{0) to second-order in the rotation rates. Starting from 
the obvious definition 



nA{r, 0%^^^^g^ =0, 
and writing the second-order surfaces as 

RA{e) = Rf + 5RA{e) + o{n^) , 

we can express the second-order correction terms as 
5nA{r, 9) 

R 



SRa{0) 



.(o)V 



(100) 
(101) 

(102) 



where primes again denote radial derivatives. We define 
the ellipticity of a surface in terms of the radii of the equa- 
tor i?cqu and of the pole i?poie as 



" R 



•cqu 



To second-order, this means that 



ba 



2R 



(0) 



Sua,: 



Aoy 



(103) 



(104) 



Here and in the following we will assume that the two 
surfaces coincide in the static case, i.e. Rn — Rp —I, 
in order to simplify the discussion. 

The respective moments of inertia Ia to second-order 
can be written as 



Ia^I 



(0) 



5Ia + 0{[f 



(105) 



where the second order correction term can be seen to be 



(106) 



(97) SIa == m / m^SuAir, 9) dV , 



which reduces to the explicit integral 



5Ia = 



Sttto 



5nAfl{r) - -SnA,2{r) I r 



^dr. 



(107) 



3 JO 

The Kepler rotation rate J^k is reached when one of the 
two fluids spins at the rate of a test particle in Keplerian 
orbit around the equator of the star, i.e. at 9 = tt/2 and 
f = -Ra(71'/2), where we assume that A represents the 
"outer" fluid at the equator. This means that 

{dMr,n/2)~rn'^)^^^^^^^^^^Q. (108) 

Due to the rotation-induced change in the shape of the 
star, the Kepler rate also has a 0{Q^) correction, and so 
we write 

^kA = ^lo) + s^kA + o{n''), (109) 

where r2(o) is the Kepler rate around a spherical star, i.e. 

4 

3' 

Using the second-order expansions for Ra{9) and $ this 
becomes 



V^l^ ^ $(o)'(l) = i;^Gp 



(110) 



5VL 



K,yl 



3$(°)'(1) 5Ra{^/2)+5<^'{1, ^/2)+0{n^) , (111) 
we can reformulate this 



and further using (102) and (|66| 
as 

3 



ml. 



fc^(i) 



5nA(l,7r/2) + 5$'(l,7r/2). 



(112) 



The Kepler rotation rate fix determines the maximal ro- 
tation rate of the respective fluids before mass-shedding 
will occur at the equator, i.e. if ^ denotes the "outer" fluid 
at the equator, then 

^A < nK,A ■ (113) 

Thus (|ll2| ) is a quadratic expression that determines the 
respective maximal rotation rates of the two fluids. Given 
the results shown in Fig. [^ we would expect the above 
equations to determine the Kepler limit to within 10-20%. 
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8. An analytic solution 

We have now completed the description of our slow- 
rotation formaUsm for two-fluid systems. Given any suit- 
able equation of state (which must provide all the relevant 
parameters, eg. pertaining to entrainment) equations ( |73[) 
and (JTJ) can be solved to produce a rotating configura- 
tion. For a typical equation of state (EOS) the calculation 
will obviously require numerics. Although we have written 
a code that solves this problem, we will not discuss such 
numerical results here. Instead we focus our attention on 
a somewhat surprising fact: It is possible to find an ana- 
lytic solution to our equations, including the general case 
where the two fiuids are rotating at different rates, for 
a particular class of EOS. This model EOS corresponds 
to constant structure functions Sxy and 0^ (and there- 
fore also constant coefficients fc, E^^ ...). As we will argue 
below, the resultant model is reasonably realistic and we 
expect it to prove useful in future studies of the dynamics 
of superffuid neutron stars. Before deriving the analytic 
solution, we shall analyze this particular class of EOS and 
assess its physical relevance. 

8.1. The "analytic" equation of state 

Within the present slow-rotation approximation any equa- 
tion of state can (quite generally) be expressed as 



£{nn, Tip, A^) = f("'(7i„, 7ip) + a(")(n„, 7ip) A^ 
+ O(il^), 



(114) 



since the relative velocity A is 0{^. We recall that the 
derived "structure functions" E^"^ , k etc. of ( |60| ) and 
(p3) are determined in terms of the two basic structure 
functions, the "density structure" SxY{r), cf. ([55|), and 
the "entrainment structure" l3^{r), cf. (p6[). In the slow- 
rotation expansion ( 114 ), they become 



Sxy — 



\dnxdnY ) „ 



and 






''X 



(115) 



The condition of constant structure functions determines 
the following class of equations of state: fj 



£ = -nx A 



XY 



Uy 



A^B^ 



nx ; 



(116) 



with constant coefficients A^^ (symmetric) and B-^ . 
Therefore (IllSl) leads to the identifications 



S: 



'XY 



(A-) 



XY 



and P^ ^B^ 



(117) 



This class of EOS is the natural generalization to 
the two-fiuid case of the "analytic polytrope" for sin- 
gle fluids {£ oc p^) for which the rotational corrections 
can be found analytically (Chandrasekhar and Lebovitz 



It is possible to add linear terms to £^ ' and a constant 
to q'"', but this does neither change the resulting structure 
functions nor the static solution, and such terms have therefore 
been omitted. 



1962). This model EOS also contains the simple case of 
the sum of two such polytropes without entrainment, i.e. 



£ ^A" 



APPrip, for which the corresponding analyti- 



cal slow-rotation solution has already been found by Prix 
(|199§). 

The major novelty of (|116| ) concerns the inclusion of 
entrainment, and it is obviously important to evaluate 
how well this model approximates "realistic entrainment." 
From the relation (n9h between entrainment and effective 



masses we see that the equation of state (116) is simply 
characterized by a constant effective mass -mP* . We noted 
earlier that nuclear physics calculations of the effective 
proton mass imply that this is not an unreasonable ap- 
proximation, even in the extreme context of neutron stars. 



This suggests that, not only does (116) have the attrac- 
tive feature of leading to an analytical solution for slowly 
rotating conflgurations, it also appears to provide a rea- 
sonable approximation which should correctly reflect the 
main qualitative features of neutron star matter (in par- 
ticular the entrainment). 



Another novel feature of (116) is the off-diagonal term 
A"P, which leads to the term 5np. It will turn out that the 
combination 



Snp dfin/dn 
cr = — 



5pp d^n/dUn 



(118) 



plays an especially important role in determining the 
structure of our neutron star model. For a reasonable EOS, 
the matrix Sxy has to be positive definite, which implies 
the constraint 



t^np "^ <-*nn<-*pp • 



(119) 



This then determines the range of acceptable values for 
a. For realistic situations the proton fraction will be 
Xp < 0.5, in which case one finds Spp < Snn- Therefore we 
can be sure that in this case ( 119 ) is satisfied if cr € [—1,1]. 
We have checked the reasonableness of this range by de- 
termining a for a more realistic EOS due to Prakash et 
al ( 1988 ). For the range of densities and proton fractions 
to be considered here, we find that a typically decreases 
monotonically from values slightly larger than one at nor- 
mal nuclear matter density Pnuc, to values around —0.4 
at about lOpnuc, with the zero occurring at around 2pnuc- 
Thus, it seems that restricting a constant a to the range 
(T e [—1,1] is reasonable. Finally, we note that a is inti- 
mately related to the so-called "symmetry energy" term 
in the Prakash et al ( |1988[ ) EOS, which is designed to 
vanish whenever there are equal numbers of neutrons and 
protons. 

8.2. The static solution 

We begin by discussing the static configuration, which is 
the solution to (^) and (^) . The chemical potentials ^ 
for the "analytic" equation of state (116) are given by 



// 



(0) = A^^n(°) 



Y 1 



(120) 
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which vanish at the surface of the star, and therefore the 8.3.1. The 
static constants of integration C'^(o) are determined from 
(Eq) at the surface r = 1, i.e. 

Eq. (MS) can now be expressed in the form 



The solution for the I — 2 component is independent of the 
chosen stellar sequence and is determined by the boundary 
(121) condition ( pi|) alone. This yields 



$f^M 



E- 



ixy 



4")(r) = fcx ^^•'^(l) 



$iu)(r) 



(122) 



5 Jb/2{rTT) 
71 \/^ 



(129) 



Inserting this into the I = 2 component of (p2) gives the 
and therefore the ratio of the respective densities in the corresponding total mass density coefficient 
static configuration with this equation of state is seen to 



be a constant, i.e. 



pf^M 



-E- 



XY 



nf{r) 



5 J5/2irTT) 

72 V^ 



(130) 



n^^\r) 



(123) while the respective particle density coefficients are found 
from (pll) 



This implies that the two fluids share a common outer sur- 
face at r = 1. Inserting the densities (122) into the second 
static equation (|46| ) results in the Lane-Emden equation 
for the total mass density p'"^ : 



XY 
^A.2 



kAE 



XY 



5 ^5/2 (^tt) 
V2 V^ 



Efjr\ 



(131) 



which has the well known solution 



P^°^M = po 



sm 



{r^/k) 



"\jk 



where po = 1 



(124) 



(125) 



8.3.2. I = 0: FCD sequence 

The Fixed Central Density solution for the I — compo- 
nent is determined by ( |9l| ) together with (pS) . In this case 
one finds that 



^^(r) = 



-B^^ /3n/2 Ji/2(r7r) 



The requirement that the density vanishes at the surface 

r — R (where here R — I) determines the constant fc, 

which is related via ( |63| ) to the symmetric structure matrix From the I = 

SxY, namely 



4-3 



k — SxY'rn m 



(126) 



component of (^ 



I) we have 
-6 



(132) 



(133) 



Integrating the density (|125|), we obtain the following re- and from (|l]) it follows that 
lation between the total mass M, radius R and the central 
density po for this EOS, 






M == -po^ 



(127) 



kAE^"^ h>^/2 Ji/2{rTi) 



E'^Jr^ 



(134) 



8.3. Analytic solution for the slow-rotation coefficients 



Given the "analytic equation of state" (|116D , the struc- 
ture functions are constant, and so are the coefficients in 



8.3.3. / = 0: FM sequence 

The Fixed Mass solution for the I = component is de- 



luie luiii^uiuiis aie uuiisuaiiu, cuiu su aie uiie uueiin;ieiius in /FTTn /PTT 

the differential equations (|| and (@). Therefore one can termined by {^ together with (^. One finds 



write their general regular solutions as 
Ji/2{rVk) , E^^ 



aXY 

^0 



tti^SabC^ 



IT 
.XY 



<^,^{r) 



E 



XY 



^ 



Ji/2{r'^) 



1 



(135) 



+ 



and from the / = component of (^) it follows that 



cl>f^(r) = A 



E 



XY J5/2{ry/k) 

\/r k 



XY 



(128) 



XY 



where Ji/2{x) and J^/2{x) are the standard Bessel func- 
tions. From now on we will set fc = tt^ in accordance "-2,^0 
with the constraint (126). The constants ^^^, ^2"^j ^'^d 
the C^'^^ are determined by the boundary conditions dis- 
cussed in section O. 



while (pll) leads to 

kAE""^ 



.W"t"^ 6 



(136) 



+ EV ( r' 



V2. 

3 

5 



Ji/2{r-K) 



r^- 



(137) 
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8.4. Physical parameters of the two- flu id star 



The analytic solution (129)-( |l37| ) has been expressed en- 
tirely in terms of the matrices i?;^^ and E^^ , and we will 
now discuss the explicit form of these coefficients in terms 
of the physical parameters describing the configuration of 
the two-fluid star. 

The symmetric matrix Sxy has three degrees of free- 
dom. One of the se is already determined by the static 
radius constraint (126) and another is associated with the 
a parameter discussed earlier. A further constraint comes 
from the proton fraction Xp, which we define as 



,(0) 



,(0) 



,(0) 



e [0, 1] 



(138) 



and which is constant throughout the static star due to 
(123). In order to simplify expressions, we take the proton 
and neutron masses to be approximately equal, i.e. 

m" « mP « m . (139) 



Using ( |126| ) and ( |123| ), the kx are expressible as 



fen = — (1 



and 



rCp — Xp 

m 



(140) 



The "density structure" Sxy can be expressed explic- 
itly in terms of these parameters as 



JXY 



{(1 - Xp) + ct(1 - 2a;p)} XpCr 



.(141) 



771^(1 + a) \ -^p^ -^p 

From this expression we see that, even though a —> —1 
seemed to be allowed from the "thermodynamical" point 
of view, it leads to a singular behaviour in the slow- 
rotation expans ion. O f course, given the fact that the 
Prakash et al ( 1988 ) equation of state indicated that 
a > —0.4 for all reasonable neutron star densities, cf. 
Sec. 8.1, we do not expect this singularity to have any 
physical significance. 

As we have seen earlier, cf. (|l9|), the entrainment a 
can be completely described by the effective proton mass 
mP*, which for our EOS ( 116 ) is a constant, and so we 
can choose without loss of generality 



/3" = , and 2/3P = me , 



(142) 



where we used dimensionless entrainment coefficient e de- 



fined earlier in (20). 

With the "structure functions" S^^ and /3^ expressed 
completely in terms of the proton fraction Xp, the "symme- 
try energy" dependence a and the entrainment coefficient 
£, the explicit expressions for E^ 

^2 



ixy 



are 



^XY 



^n 


3m(l + cr) 




/ {1 — Xp -1- o-(l — 2a;p — X 




\ Xpae 


and 






TT^Xp / (cr - e) e 
3TO(l + cr) ^ e (l-£ 



p£)} 



Xpcre 
a;pcr(l ~e) 



143) 



(144) 



while the coefficient E of (p3) is found to be 



^XY ^J^ f ^-Xp{l + S) XpE 

3 \ XpE Xp(l-£) 



(145) 



It is interesting to note that in the case of co-rotation the 
corresponding terms in the analytic solution become 



ilxEj^ Hy — — ^^ 



and 

nxE>'^nY = ^n\ 



(146) 



(147) 



and are therefore seen to be independent of not only the 
entrainment e, but also of a. 

8.5. Explicit results for ellipticity, moments of inertia, 
and Kepler rotation 



Using the analytic solution of Sec. 8.3 with the physical 
parameters of Sec. 8.4, we obtain the following explicit 
expression for the ellipticities (104) 



e.= 2(l 



15 

72 



^ I rtxE'^'^nY 



?(0) 



37r2 

2kA' 



^^xEV^y 



where we have set R\ = 1. 



For the moments of inertia, integrating (|10 
analytic solution leads to the explicit result 



5Ia 

A") 



aQxE^^^Y 



3b 

kA 



HxEa iiy 



with coefficients 



TT 

175 



b = 



3^fa 
175(7r2 - 6) 



(148) 



for the 



(149) 



(150) 



(151) 



At the formal level, this result is equivalent to that found 



in earlier work without entrainment (Prix 1999). However 



the two results differ at a deeper level since the actual ma- 
trix elements of E^^ and E^^ contain information about 
entrainment and the "symmetry energy" contribution a. 
For the analytical solution we can express the correc- 



tion to the Kepler rotation rate (111) as 



^f^2 6 /I 3 \ ^ j^xY^ 27 ^ ^^y 



'K,yl 



lOkA 



nxEi'^YA'^^'^) 



where A denotes the "outer" fiuid at the equator. 

If we take the simple case of the two fluids co-rotating 
at the maximal rotation rate, i.e. $!„ = ^p = ^k, this sim- 
ply leads to 



nK = 



n 



(0) 



0.69 r2(o) w 0.80 VttG^, 



(153) 
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independent of the parameters of this EOS. This is 
about 20%-25% too high compared to the numerical re- 
sult r^K ~ 0.64:y/TrGp. This discrepancy originates in the 
underestimate of the equatorial radius when approach- 
ing SIk, which is illustrated in Fig. P. The "failure" of 
the slow-rotation expansion to accurately determine the 
Kepler rotation limit has already been pointed out by, for 



example, Salgado et al. (1994) 



9. Exploring the roles of relative rotation, 
entrainment and "symmetry energy" 

In this section we will use the analytic solution to the 
superfluid slow-rotation problem to explore the effects of 
relative rotation, entrainment, and "symmetry energy" on 
the distribution of matter, the Kepler limit, ellipticity, 
and moments of inertia for a fixed mass star. The re- 
sults we present are for a particular stellar model with 
mass IAMq and radius 10 km. The proton fraction is 
taken to be Xp = 0.1 and we choose the neutron rota- 
tion rate to be (except in Figs. Q and 0) f^n = 0.1. This 
value corresponds to (in our units) the rotation rate of 
the fastest known pulsar (1.6 ms). The three parameters 
that will be varied are a, e, and the relative rotation rate 
Qp/rin. For the first two parameters we will consider only 
the values a = -0.5,0,0.5 and e = 0,0.4,0.7, while the 
relative rotation rate will be assumed to lie in the range 
—2 < rip/rin < 2 (which allows the protons to be counter- 
rotating with respect to the neutrons) . The values chosen 
for a are in accordance with the discussion in Sec. 8.1. We 



also recall that the best current estimates for entrainment 
imply a range of 0.4 < e < 0.7. This means that our cho- 
sen models correspond to the expected upper and lower 
limits, as well as the no-entrainmcnt limit (which provides 
a useful reference). 

9.1. Effects on the "local" structure 

First we focus on the "local" structure by examining how 
the particles get redistributed throughout the star because 
of rotation. In Fig. g we show plots of the rotationally- 
induced corrections to the neutron number densities in the 
equatorial plane and along the polar axis. We show plots 
for (T = only because other values lead to similar results. 
An inset is provided to show that the effect of entrainment 
is small but not completely negligible. Recall that we are 
considering the fixed mass solution, so the density at the 
center of the star will decrease as the star is spun up. 
This also explains why the density correction is negative 
along the polar axis. Furthermore, the density corrections 
in the equatorial plane are positive near the stellar surface 
since rotation causes the matter to bulge. Fig. H shows the 
corresponding proton particle number density corrections. 
Unlike in the case for the neutrons, we have added a graph 
for (J = 0.5. From Fig. y we notice that the effect of chang- 
ing a and e is more pronounced for the protons. That this 
should be the case is easy to understand since the proton 
fraction is small and only 10% of the mass of the star is 



^ 2 



•3.9 
■4.1 
•4.3 






e = Ji/2. 



e = 
e = 0.4 
£ = 0.7 




-2 



Fig. 2. Plots of the radial profiles of the neutron density 
corrections, normalized by the static number density, in 
the equatorial plane {6 — 7r/2) and along the polar axis 
(61 = 0) for cr = 0, e = 0,0.4,0.7, and f7„ = 0.1 and 
rip ~ 0.25r2n- In this, and the following figures, the dotted 
lines correspond to £ = 0.7, dashed have e — 0.4, whereas 
the solid lines are for e = 0. 



in the protons. The relative effect of a changing e is more 
apparent in the proton plots, then, because the absolute 
magnitude of the proton number density cor recti on s are 
alwa ys 10 % of those of the neutrons. From ( |l37|) , ( |l43| ) 
and ( 144 ) we see that the main difference between the 



two cases is that while changes in e and a affect E^"^ 
at leading order {^ Xp), they only lead to higher order 
corre ctions to E^^ (since there are also terms of 0(1) in 
(143)). We can also see from the right-hand panel of Fig. |^ 
that changing a from to 0.5 results in substantial modi- 
fications. This shows clearly that the "symmetry energy" 
can play an important role in determining the rotational 
configuration of a two-fluid star. 

Fig. y provides a different perspective on the 
rotationally-induced redistribution of the particles. Here 
we show isodensity surfaces in the (r, 0)-plane, for e = 0, 
(7 = (in the left panel) and e = 0.7, a = —0.5 (in the 
right panel). Note that we have considered slightly larger 
values of the rotation rates (just below the Kepler limit) in 
order to exaggerate certain effects. For a given density, we 
compare the isodensity curves for a non-rotating star to 
the neutron and proton isodensity curves for the rotating 
model. In both panels we see that the neutron and proton 
curves actually intersect each other near the surface of the 
star. That is, along the equator the neutrons actually ex- 
tend further than the protons, while the protons extend 
further along the rotation axis. This means that, the ro- 
tational configuration of the protons is, in fact, prolate. 
We note that this effect has been exaggerated in the pan- 
els, because of the higher rotation rates, but it happens 
also for lower rotation rates (as used in the earlier figures). 
Near the center of each panel, we see that the neutron and 
proton surfaces no longer intersect. This is no doubt due 
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Fig. 3. Plots of the radial profiles of the proton density corrections, normalized by the static number density, in 
the equatorial plane (Q = it/2) and along the polar axis {9 = 0) for a = 0,0.5, e = 0,0.4,0.7, and fin = 0.1 and 

Vtp = 0.25f]n. 
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Fig. 4. Neutron (thick line), proton (thin line) and static configuration (dashed line) isodensity curves in the meridional 
plane, for cr = 0, -0.5, e = 0,0.7, and 17n = 0.15 and Vl^ = 0.25^1- 



to the fact that the centrifugal forces are smaller closer to 
the center of the star. 



9.2. Effects on the "global" structure 

Now we focus our attention on the roles of a, e, and the 
relative rotation rate rip/rin in determining macroscopic 
properties of the star; in particular, the proton and neu- 
tron ellipticities, their respective moments of inertia, and 
the Kepler limit. As before, we will keep the mass of the 
star fixed. 

Fig. ra illustrates the neutron and proton ellipticities 
as functions of the relative rotation rate, and Fig. O gives 
their moments of inertia (normalized by their static val- 
ues). There is an obvious quadratic behavior in each plot 
due simply to the slow-rotation expansion. As well, the 
intersection of all curves at flp/iln — 1 occurs because 
the protons then co-rotate with the neutrons and the sys- 



tem is behaving as a single fluid. Notice that entrainment 
has the largest influence when the neutrons and protons 
counter-rotate. This is easily understood as a consequence 
of the basic fact that the entrainment parameter repre- 
sents the way that the equation of state (as represented 
by the energy functional) depends on |t)n — Vpp. The pro- 
tons are in general much more affected by changes in the 
various parameters than the neutrons, again due to the 
fact that the neutrons carry the bulk of the mass of the 
star. Perhaps most interesting is the effect of both the en- 
trainment parameter e and "symmetry energy" parameter 
(7 in determining the minima of the curves. For both the 
protons and the neutrons, wc sec that an increase in e for 
a fixed value of a leads to a deeper value for the minimum. 
Decreasing the value of a causes the minima to become 
even deeper. In particular, we see from the left-most panel 
in Fig. B that the minima have become deep enough that 
the protons can be prolate (i.e. have negative ellipticity) 



16 



R. Prix et al.: Slowly rotating superfluid Newtonian neutron star model 



even though they rotate in the same direction as the neu- 
trons. Finally, we note that as a is decreased, the neutron 
ellipticities go from having minima in the right-most and 
center panels, to having maxima in the left-most panel. 
That is, if the absolute value of the relative rotation could 
be made large enough the neutron fluid could also become 
prolate. 

It is interesting to compare the curves for the elliptic- 
ities in Fig. and the corresponding moments of inertia 
in Fig. pi for example for cr = 0, e = 0.4 (middle panel, 
dashed line) for the protons. For fip = we observe that 
the proton fluid is nearly spherical (ep « 0), and still its 
moment of inertia is higher than of the static (spherical) 
configuration 51-p > 0, and this is simply because the mass 
distribution has been shifted further away from the ro- 
tation axis even though the surface itself is (nearly) un- 
changed in this case. This can also be clearly seen in the 
density correction in Fig. y, for example for a = 0, e = 0.7 
and rjp/rin = 0.25 (left panel, dotted line), which nearly 
vanishes at the surface r = 1, and is to be compared to 
the corresponding global quantities in Fig. ^ and Fig. ^j. 

Finally, results for the Kepler, or mass-shedding, limit 
are shown in Fig. M. To understand these results one must 
appreciate that there is a subtle difference with the single- 
fluid case: In our case the two fluids can rotate indepen- 
dently at different rates. Thus, one of the fluids typically 
extends beyond the other, in particular at the equator. 
Since the Kepler limit is deflned by the outermost fluid 
at the equator, we can use Eq. ( |109| ) in the following 
way: When the neutrons are outermost, set ^ = n and 
fin = ^K.n and solve the resulting quadratic for ^x.-a 
as a function of the ratio fin/fip, and vice versa in the 
case when the protons extend beyond the neutrons. In 
Fig. we show the resultant solutions over the entire 
range of the relative rotation rate. The Kepler rate is easy 
to determine, however, because it is given by the neu- 
tron curves when fin/fip > 1, and the proton curves when 
fin/fip < 1. Of course, the various curves always intersect 
when J7n/fip = 1, the case that corresponds to corotation 
of the two fluids. For the case of cr = e = 0, we find re- 
sults in good qualitative agreement with the relativistic 
study of Andersson & Comer ( 2001a ). In particular, we 
see that the Kepler limit changes lit tle wh en fin > fip. As 
explained by Andersson & Comer ( 2001a ), this is due to 
the fact that the neutrons make up 90% of the mass of the 
star, and the star is behaving like a single-fiuid star with a 
small proton component. When fip > fin the Kepler limit 
increases as fin is decreased. Again, this is natural because 
the neutrons still dominate the mass of the star, and the 
Kepler rate is simply approaching the non-rotating star 
limit. 



10. Conclusions 

We have developed a formalism for modeling slowly- 
rotating Newtonian superfiuid neutron stars incorporat- 
ing entrainment. We have used a two-fluid description, 
where one fluid is the superfluid neutrons and the other 



is a charge-neutral conglomerate of the remaining con- 
stituents. A detailed discussion of the relation between 
entrainment and nuclear physics calculations (i.e. equa- 
tions of state) was given. Using an equation of state that 
is quadratic in both the mass-densities and relative veloc- 
ities of the fluids, we found that an analytic solution to 
the slow-rotation equations could be obtained. This solu- 
tion is the natural extension to the two-fluid case of the 
E (x p^ polytrope in the single fluid case (which has proven 
to be very useful for understanding the properties of or- 
dinary fluid neutron stars). We used the analytic solution 
to explore effects due to relative rotation, entrainment, 
and "symmetry energy" on both the "local" and "global" 
properties of a fixed-mass star. An unexpected result was 
that the "symmetry energy" parameter had as much im- 
pact on the rotational equilibria as the entrainment pa- 
rameter. 

Our ultimate goal is to study the modes of oscillation 
of both Newtonian and general relativistic slowly rotating 
superfluid neutron stars. We believe that the formalism 
and analytic solution presented here will be valuable in 
reaching this goal. In particular, the inclusion of entrain- 
ment is absolutely necessary in determining how the domi- 
nant dissipative mechanism (the so-called mutual friction) 
of rotating superfluid neutron stars affects the gravita- 
tional radiation emitted from unstable modes. 
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Appendix A: An alternative derivation of the first 
integrals of motion 

A more elegant way of finding the first integrals of motions 
^ consists of using the identity 



v^Vj Pj = CuPi - Pj VjW-' , 



(A.l) 



in terms of the Lie derivative C^, which allows one to 
rewrite Euler's equations (|2|) in the form 

(A.2) 
(A.3) 



dtp^ + £„^pp_v.(pP + «^pP)=o. 

Stationarity {dtpf = 0), uniform rotation (wj^- = fixv') 
and axial symmetry {Cippf = 0) then directly results in 
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Fig. 5. The neutron and proton ellipticities as functions of the relative rotation Op/iln, for a = —0.5,0,0.5, e 
0,0.4,0.7, and ^i = 0.1. 
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and inserting (|9|) this explicitly becomes 
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